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We study the physics of interacting spin-1 bosons in an optical lattice using a variational 
Gutzwiller technique. We compute the mean-field ground state wave-function and discuss the evo¬ 
lution of the condensate, spin, nematic, and singlet order parameters across the superfiuid-Mott 
transition. We then extend the Gutzwiller method to derive the equations governing the dynamics 
of low energy excitations in the lattice. Linearizing these equations, we compute the excitation 
spectra in the superfiuid and Mott phases for both ferromagnetic and antiferromagnetic spin-spin 
interactions. In the superfiuid phase, we recover the known excitation spectrum obtained from Bo- 
goliubov theory. In the nematic Mott phase, we obtain gapped, quadratically dispersing particle 
and hole-like collective modes, whereas in the singlet Mott phase, we obtain a non-dispersive gapped 
mode, corresponding to the breaking of a singlet pair. For the ferromagnetic Mott insulator, the 
Gutzwiller mean-field theory only yields particle-hole like modes but no Goldstone mode associated 
with long range spin order. To overcome this limitation, we supplement the Gutzwiller theory with 
a Schwinger boson mean-field theory which captures super-exchange driven fluctuations. In addition 
to the gapped particle-hole-like modes, we obtain a gapless quadratically dispersing ferromagnetic 
spin-wave Goldstone mode. We discuss the evolution of the singlet gap, particle-hole gap, and the 
effective mass of the ferromagnetic Goldstone mode as the superfluid-Mott phase boundary is ap¬ 
proached from the insulating side. We discuss the relevance and validity of Gutzwiller mean-field 
theories to spinful systems, and potential extensions of this framework to include more exotic physics 
which appears in the presence of spin-orbit coupling or artificial gauge fields. 


I. INTRODUCTION 


Recent progress in ultra-cold atoms has made it possi¬ 
ble to study strongly correlated phenomena in a number 
of different contexts, which have no natural analog in con¬ 
densed matter physics. One such example is spinless and 
spinful many-body bosonic systems, whose rich physics 
has been experimentally explored in the continuum and 
in optical lattices M- Large spin systems offer inter¬ 
esting possibilities to study the interplay between com¬ 
peting/complimentary orders at zero and finite temper¬ 
atures such as single-particle and pair superfluidity, spin 
and liquid crystallinity, all of which can be probed using a 
variety of experimental tools Recently, attention 

has turned to the physics of bosonic and fermionic sys¬ 
tems in the pres ence of spin-orbit coupling or artificial 
gauge fields [lll - ll^ . Spin-orbit coupling and artificial 
gauge fields introduce degeneracies in the single-particle 
spectrum, which tends to frustrate the usual Bose con¬ 
densation, setting the stage for the appearance of ex¬ 
otic ordered phases even at the mean-field level, such 
as striped Bose condensates which spontaneously break 
translational symmetry [l9l - l^ or magnetized spin stripe 
phases which spontaneously break time-reversal symme¬ 
try [^ . Furthermore, single particle degeneracies can 
amplify the role of quantum fluctuations leading to chiral 
superfluids or bosonic phases with topological order, 
which resemble the integer and fractional quantum Hall 
effect [ 2 ^ of electrons in a magnetic field. 
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With these exciting experimental advances, it has be¬ 
come increasingly important to develop theoretical meth¬ 
ods which are sophisticated enough to describe the mul¬ 
titude of order parameters and broken symmetry phases 
that can potentially occur in these interesting interact¬ 
ing systems, even at the mean-field level. In addition, it 
is essential to first gain significant insight into the var¬ 
ious forms of order such systems can develop before in¬ 
troducing the next level of complexity through artificial 
gauge fields or spin orbit coupling. With this in mind, 
we perform a variational Gutzwiller study of the mean- 
field physics of a spin-1 Bose gas in an optical lattice, 
in the absence of spin-orbit coupling. We highlight the 
key virtues and limitations of this method in describing 
spinful bosonic systems, by giving a comprehensive ac¬ 
count of the static and dynamic properties of the spin-1 
Bose gas in an optical lattice. We then supplement this 
method by a Schwinger boson mean-field theory to cap¬ 
ture the gapless Goldstone mode in the ferromagnetic 
Mott phase. We discuss extensions of this approach to 
the spin-orbit coupled, large spin problem. 

The bosonic Gutzwiller mean-field technique, intro¬ 
duced by Rokhsar and Kotliar [ 2 ^, describes the mean- 
field physics of the Bose Hubbard model , which is re¬ 
alized by trapping bosons in a deep optical lattice ii- 
Known to be exact in infinite dimensions, the technique 
captures local physics by decomposing the full Bose Hub¬ 
bard Hamiltonian into a sum of on-site Hamiltonians cou¬ 
pled by mean-fields. The transition from the superfiuid to 
the Mott insulator is then obtained by self-consistently 
solving for where the mean-field order parameter van¬ 
ishes. Within the Mott lobes, the Gutzwiller Hamilto¬ 
nian is therefore purely local {i.e., identical to the zero 
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hopping limit), and hence does not a priori capture any 
correlations. In the spinless case, where exact Quan¬ 
tum Monte Carlo (QMC) simulations are possible (^ . 
the Gutzwiller method is known to correctly capture the 
qualitative features of the phase diagram. We will apply 
the Gutzwiller technique to spin-1 bosons in the current 
work. 

The equilibrium physics of this model for spinless 
bosons is well known, and recently it has been extended 
to include spinor systems, such as the one we study here 
(29l - l^ . For spinless bosons, this technique has recently 
also been extended to finite clusters, where the physics 
within an m X n plaquette is solved exactly, and the pla- 
quettes are coupled by mean-fields. Systematic studies 
using this cluster Gutzwiller method have shown that for 
relatively small cluster sizes, quantitative agreement is 
obtained for the location of the phase boundaries with 
numerically exact Quantum Monte Carlo studies [35l| . 
Importantly, such cluster methods offer a comparatively 
numerically efficient way to systematically incorporate 
correlation effects and perform non-equilibrium dynam¬ 
ics, which is highly relevant to ongoing experiments on 
strongly correlated bosons [Sbl - I^ . 

In this paper we provide a systematic study of the 
physics of the spin-1 Bose Hubbard model. We dis¬ 
cuss the equilibrium theory and highlight the advan¬ 
tages and disadvantages of the Gutzwiller approach in 
correctly capturing the spin physics in an interacting 
bosonic system. Although the equilibrium phase bound¬ 
aries for this model have been established by several au¬ 
thors m, HO, 113 , a systematic study of the various order 
parameters present in the spin-1 Bose Hubbard model 
has been lacking. As this is essential in forming a com¬ 
prehensive understanding of each phase, a new feature of 
our work is a focus on the evolution of the various order 
parameters, such as the spin, nematic and singlet order. 
As these order parameters can be directly measured in 
experiments, such a study is a necessity in the interpre¬ 
tation of the experiment and the understanding of the 
spinful Bose-Hubbard quantum phase diagram. 

Recently, several authors (dol - ld^ have extended the 
Gutzwiller technique to capture dynamics, by lineariz¬ 
ing the mean-field equations of motion about the su¬ 
perfluid and Mott insulating ground states. Indeed the 
linearized theory correctly captures the well-known Bo- 
goliubov modes in the superfluid, the gapped particle- 
hole like modes in the Mott insulator, and correctly de¬ 
scribes qualitative aspects of how these modes evolve 
across the superfluid-Mott phase boundary. Here we 
extend the time-dependent Gutzwiller framework devel¬ 
oped by Krutitsky and Navez (dj to spin-I Bose Hub¬ 
bard model, and compute the low energy spectrum across 
the entire superfluid-Mott phase diagram for ferromag¬ 
netic and anti-ferromagnetic spin-dependent interactions 
( 43 . In particular, we show that the Gutzwiller ap¬ 
proach by itself, fails to fully capture the spin physics 
in the Mott insulator, and therefore we supplement this 
theory with a Schwinger boson mean-field theory, which 


captures inter-site magnetic fluctuations. Our combined 
Gutzwiller -|- Schwinger boson approach thus more or less 
fully captures the low energy mean-field properties of the 
spin-1 Bose Hubbard model, and sets the stage for stud¬ 
ies of more complicated correlated bosonic Hamiltoni¬ 
ans, which include spin-orbit coupling or artificial gauge 
fields, which we leave for future study. 

The rest of this paper is organized as follows: in 
Sec. H, we present the Gutzwiller mean-field equations, 
and define the various order-parameters which we use to 
compute the ground state properties and the superfluid- 
insulator phase boundaries. In Sec. HI, we present the 
Gutzwiller phase diagrams for ferromagnetic and anti¬ 
ferromagnetic interactions focussing on the evolution of 
the spin, nematic and singlet order parameters across 
the entire phase diagram, which are absent in the spin¬ 
less case. In Sec. IV, we linearize around the Gutzwiller 
ground state to produce the equations for the low energy 
collective modes. In Sec. V, we present the excitation 
spectra for ferromagnetic and anti-ferromagnetic inter¬ 
actions in the superfluid and Mott phases. In Sec. VI, 
we discuss the relevance of our results to experiments and 
in Sec. VH, we provide a summary of our results and dis¬ 
cuss directions for future work. Readers uninterested in 
the detailed theoretical derivations of the equations may 
skip Sections H and IV, which are technical in nature. 


II. GUTZWILLER MEAN-FIELD THEORY: 

STATICS 

In this section we outline the Gutzwiller mean-field 
theory which we use to compute the ground state energy, 
wave-function and order parameters in the following sec¬ 
tions. We remark that the Gutzwiller mean-field equa¬ 
tions have been obtained previously by several authors 
[13, HO, 113 and are only reproduced here for complete¬ 
ness (and for providing a background, as well as a context 
for the new results obtained by us for the spin-I bosons). 
Our starting point is the generalized Hamiltonian of the 
spin-1 Bose gas in an optical lattice [13 

H - pN = -t ^ (al^aja + H.c^ +^'^ni{hi- 1) 

i i 

where {i,j) denotes nearest neighbor sites, a chemi¬ 
cal potential p, and we have defined fii = 

S = p apaip , where the T are a vector of spin-1 
matrices. While the mean field equations are generic for 
any dimensionality and lattice geometry, in the following 
we focus on a three dimensional cubic lattice. 

To the derive the Gutzwiller equations, we proceed by 
writing the tight binding model at the mean field level 
by treating each neighbor of site i using the mean field 
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the Hamiltonian (Eq. [T]) with respect to the mean-field 
ansatz (Eq. |3|). The terms diagonal in particle number 
(2) yield 


Following the standard Gutzwiller mean-field theory, we 
write the ground state wave function in the Fock basis as 
a direct product over single site wave functions as 

i^-Gs) = oilr !</>*) (3) 

= E (4) 

m_i ,mo ,mi 
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where m_i,mo,mi denote Fock states in the mz = 

— 1,0,1 state respectively, and we have assumed the The spin operator yields off diagonal terms in Fock 

^m_imomi to be site independent. We will now evaluate space. In particular we obtain 


(<(.,1^(82-2h.)|<^.) = C/2 ^ 


m_i 




-m_i] -I-ruimo-b (5) 


^m_imomi^(m_i-l)(mo-|-2)(mi-l)'\/wi("lO + l)(wo -f 2 )to_i 

\m.i+i)(mo-2)(mi+i)V{'rni + l)mo(mo - l)(m_i -b 1) 
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Now that we have all of the two particle terms we can 
move on to calculating the expectation value of hopping 
{i.e. t) dependent terms. For the creation operator this 
yields 

^ ( 6 ) 

-^T7i_imo(mi —1) “1“ ^m_i (mo —l)?7ii 

+A( 

where z denotes the number of nearest neighbors, and for 
the destruction operator this yields 

(^il ~ t 'y = —‘2.Zt ^ ^ ^m_imomi C^) 

06 m_i,mQ,mi 

1 4“ 1 T ^m_i(mo + l)mi T 1 

Tblyn_imo (mi-t-1) T"^ . 

We self consistently determine the values of (ajp), 
by viewing the ground state expectation value of the 
Hamiltonian as A.TL.A, where .A is a vector of all 
the Am-imomi, and we just have to diagonalize H = 
{ni,no,n_i\H\m-i,mo,mi} iteratively. This is equiva¬ 
lent to determining the Am_imomi coefficients variation¬ 


ally through solving 

-('I'Gsl-fC - = 0, (8) 

and similarly for Am_imomi- At each step of the numer¬ 
ical iteration, we self-consistently determine the mean- 
field {ata), which is identical on every site in the homo¬ 
geneous system we study. In what follows, we therefore 
occasionally drop the site index when unnecessary. 

We note that the mean-field theory described here can 
be readily generalized to capture translational symmetry 
breaking phases by allowing the mean-field to vary on 
every site. However, absent spin-orbit coupling or artifi¬ 
cial gauge fields, we do not expect translational symme¬ 
try breaking phases to occur, and therefore restrict our 
study to spatially homogeneous mean fields. 


A. Order parameters 

The spin-1 system generally possesses four order pa¬ 
rameters namely, the condensate fraction, the spin, 
the nematic director, and the singlet order parame¬ 
ter, which we determine using the ground state wave- 
function computed above. The condensate fraction in 
the three hyperfine spin states can be expressed as (a) = 
















4 



t/t/o 

FIG. 1: (Color Online) Density plot showing the evolution 
of the superfluid order parameter in the ferromagnetic {U 2 ~ 
—O.lC/o) spin-1 Bose gas. The superfluid order parameter 
goes to zero at the superfluid-Mott phase boundary as shown 
in the inset. The total spin <S/n = 1 throughout the phase 
diagram. The numbers label the density in the Mott lobes. 
All transitions here are second order. 


{(a_i), (ao), (ai)} and is given simply by 

—l)mi —J*^ ■ (^) 

The spin is a single-particle vector order parameter 
given by (S) in the mean-field ground state. The ne¬ 
matic director is a two-particle tensor order parameter 
which is given by the matrix Afap = ^{SaSp -I- S^Sa), 
where a,/3 € {x,y, z} correspond to the spin-operators 
representing the x, y, and z, components of the spin. 
Diagonalizing the nematic tensor yields three eigenval¬ 
ues, the largest of which (denoted by X^f) corresponds to 
the degree of nematicity. The singlet order parameter is 
also a two-particle order parameter which measures the 
number of singlets on a given site. It is a scalar object 
obtained by computing the number of singlets in 

the mean-field ground state, where 0|^ = 

(see Refs. is the singlet creation operator on site 

i. Once the ground state wave-function is obtained, each 
of these order parameters can be readily computed nu¬ 
merically. 


III. PHASE DIAGRAM 


Bogoliubov theory of the spin-1 Bose gas [^, @ and the 
strong coupling expansion ]T^ . . We stress here that 

while our results for the phase boundaries are not new, 
and have been discussed by several authors [29l - l3^ . we 
focus on the evolution of the spin, nematic, and singlet 
order parameters across the superfluid-Mott transition, 
which has not yet been discussed in the literature and 
is essential to understand and characterize the nature of 
each phase as manifested in experiments. Importantly, 
these order parameters distinguish the spin-1 gas from 
its well studied spinless counterpart and are generically 
present in higher spin systems, such as spin-3 Chromium 
atoms 0. Understanding how these evolve in the spin-1 
Bose Hubbard model is therefore crucial to developing 
theories of larger spin systems, which are currently being 
explored experimentally |T7H4^ . 


A. Ferromagnetic interactions 

We start by discussing the conceptually simple fer¬ 
romagnetic case for U 2 < 0, which occurs for ®^Rb, 
where the system has only two order parameters, the 
superfluid order parameter and the total spin S = 

In Fig.[Tl we plot the theoretically calculated evolution 
of the superfluid order parameter, which clearly shows 
the superfluid-Mott phase boundary, where the super¬ 
fluid order parameter vanishes. The second order na¬ 
ture of the phase transition [ 2 ^, IsO, 113 is clear from the 
smooth manner in which the order parameter goes to 
zero. This is to be expected, as on the ferromagnetic 
side, the total spin simply locks to the density, and the 
system resembles a spinless gas, which is known to have 
a second order superfluid-Mott insulator transition. The 
total spin S/n normalized by the total density is equal 
to unity throughout the phase diagram. 

A curious feature of this mean-field theory is that the 
dependence on U 2 completely drops out for the n = 1 
Mott lobe, as is evident from Eq. [S] This is because, ab¬ 
sent the mean-field term, number fluctuations are com¬ 
pletely frozen out and the n = 1 Mott lobe is obtained 
by setting either mi, mo or m_i = 1, while the others 
are zero. The physics of the n = I Mott insulator thus 
has to be inferred by continuity from the superfluid side. 
For U 2 < 0, where the superfluid is ferromagnetic, the 
Mott insulator is also ferromagnetic, whereas for U 2 > 0, 
where the superfluid is polar, the Mott insulator is ne¬ 
matic. A better treatment of the ferromagnetic Mott 
insulator, which captures spin fluctuations can be done 
using a Schwinger boson mean-field theory, described in 
detail in Sec. VB.I. 


In this section we present the phase diagram of the 
spin-1 Bose Hubbard model and highlight the key virtues 
and shortcomings of the Gutzwiller approach with re¬ 
spect to other approaches, namely the weak coupling 


B. Anti-ferromagnetic interactions 

We now turn to the more interesting case of anti¬ 
ferromagnetic interactions, which are present in ^^Na. 
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FIG. 2: (Color Online) Density plot showing the evolution of the nematic (left) and singlet (right) order parameters in the 
anti-ferromagnetic {U 2 = O.lDo) spin-1 Bose gas. The nematic order parameter is computed by taking the maximum eigenvalue 
of the nematic tensor. In the singlet phases, labelled 1, nematic order vanishes, and the superfluid-Mott transition is first order, 
as indicated by the finite jump in the nematic order, as shown in the inset. In contrast, the transition into the lobe labelled II 
is a continuous second order quantum phase transition. On the right, the singlet number takes on an integer value equal to 1(2) 
in the n = 2(4) Mott lobes, corresponding to the number of local singlets per site. The inset shows the singlet order parameter 
for two different cuts across the phase diagram: the cuts traverse the n = 2 (dashed) and n = 1 (solid) Mott-superfluid phase 
boundaries. At the n = 2 Mott-superfluid phase boundary, the singlet order parameter displays a finite jump, which is evidence 
of the first order nature of the transition. 


Here the system is described by three order parameters, 
the complex superfluid order parameter, the tensor ne¬ 
matic order parameter Af /n^, and the scalar singlet or¬ 
der parameter, which measures the number of singlets 
created on a given site. 

In Fig. [21 we present the phase diagram for anti¬ 
ferromagnetic spin-dependent interactions U 2 > 0 in 
terms of the nematic and the singlet order parameters. 
To display the nematic order parameter, we compute the 
largest eigenvalue Xj\f of the nematic tensor Mjr? for 
every value of t/C/o and /x/C/q throughout the phase dia¬ 
gram. The inset shows a cut of the nematic order param¬ 
eter at fixed /i through the n = 2 Mott lobe. The labels I 
and II indicate the nature of the superfluid-Mott insula¬ 
tor transition as being first and second order respectively. 

In the superfluid phase, we recover the well known con¬ 
tinuum result, namely that the largest eigenvalue of the 
nematic tensor is Ajy —>■ 0.5. In the simerfluid, the spin, 
density and nematic orders compete [8|, which enforces 
a constraint on the nematic tensor such that the sum of 
the eigenvalues equals 1. In the even integer Mott lobes, 
the system enters the singlet phase, which is character¬ 
ized by zero total spin (S) = 0, and zero nematicity. In 
other words, all three eigenvalues of the nematic tensor 
are zero in this phase. Unlike the ferromagnetic gas, sev¬ 
eral authors have argued that the transition from the 
singlet Mott insulator to the nematic superfluid is first 
order [sol [ 3 ^ , and this has been confirmed numeri¬ 
cally in 1 and 2D through exact Quantum Monte Carlo 
simulations (Sll. l33|. 

The first order nature of the transition should also be 


evident in the evolution of the nematic, condensate and 
singlet order parameter across the phase boundary. In 
the inset on the left panel, we plot a cut through the 
nematic order parameter across the n = 2 Mott lobe, 
which clearly shows a discrete jump at the superfluid- 
singiet Mott transition. As the nematic order can be 
probed exp erimentally using optical birefringence tech¬ 
niques |50| . this is an example of a first order quantum 
phase transition that can be readily studied in the labo¬ 
ratory. Additionally, the superfluid order parameter also 
shows a discrete jump at this transition to zero, unlike 
for a second order transition. Near the odd integer Mott 
lobes, the transition is once again second order. 

Unlike the nematic order parameter, which is identi¬ 
cally zero in the singlet Mott insulator, the singlet order 
parameter (0l0) takes on an integer value equal to half 
the total particle number in the even Mott lobes. The in¬ 
set on the right panel shows two horizontal cuts through 
the phase diagram across the superfluid - n Mott tran¬ 
sition where n = 1 (solid) and n = 2 (dashed). The 
singlet order parameter is zero in the n = 1 Mott lobe, 
as expected, and precisely 1 in the n = 2 Mott lobe. Fur¬ 
thermore, the singlet order parameter shows a small but 
finite jump at the Mott-superfluid transition, once again 
revealing the first order nature of the transition. 

A first order superfluid-Mott transition usually implies 
a small but finite coexistence region where a metastable 
superfluid phase can occur in addition to a Mott insu¬ 
lator with local singlets. Such a coexistence region has 
been discussed within the mean-field context [s^ and val¬ 
idated through exact Quantum Monte Carlo simulations 
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in ID (sTj. Nonetheless, the first order transition can 
be readily probed by studying the evolution of the su¬ 
perfluid, nematic or singlet order parameters. Precisely 
figuring out whether this transition is indeed first order or 
an artifact of the mean-held approximation will demand 
more sophisticated numerical simulations (e.q. QMC) in 
two and three dimensions, which our work should moti¬ 
vate. 

For very small U 2 /U 0 , Imambekov et al. 0 hnd an 
additional hrst order phase transition within the Mott 
lobe, which corresponds to a transition from a nematic 
Mott phase to a singlet phase. This transition is not 
captured by the Gutzwiller mean-held theory, because it 
does not include any spin huctuations in the Mott insu¬ 
lator. Nonetheless, the Gutzwiller theory gets the cor¬ 
rect local spin physics at zero hopping, and extrapolates 
this wave-function throughout the entire Mott lobe. Note 
that within this mean held theory, the singlet Mott lobes 
are larger than the nematic lobes. This is because the 
tendency to form local singlets which is favored by the 
repulsive spin-dependent interactions tends to destroy su- 
perhuid order more easily, thus enhancing the Mott re¬ 
gion of the phase diagram. 


IV. GUTZWILLER MEAN-FIELD THEORY: 

DYNAMICS 

Having formed a more comprehensive understanding 
of the static ground state quantum phase diagram using 
the Gutzwiller technique, we now turn our attention to 
the low lying excitation spectrum in the spin-1 gas. Some 
of the results derived here for U 2 > 0 were recently ob¬ 
tained in Ref. Here we present some more details on 
the excitation spectrum for [/2 > 0, and also describe the 
excitations on the ferromagnetic side U 2 <0. We start by 
presenting a complete derivation of the equations of mo¬ 
tion, which describe the full mean-field dynamics about 
the Gutzwiller ground state, generalizing the earlier work 
in Ref. [4l| for the spinless Bose Hubbard model. 

To begin, we generalize the Gutzwiller wave-function 
to include dynamics 


the time dependent Schrodinger equation with respect to 
have 




+ = 0 . ( 11 ) 


Taking the variational derivative, we arrive at 

+ S+- —l)(mo-|-2)(mi— 1 ) (12) 


+ Srr^ ^ 


(m_i + l)(mo—2)(mi + l) 


— t 




T — T y/^l^rri-i tuq (mq — 1) 

- t yjm-i -I- lA 

(m_i+l)mQmi -I- ^mo + lA 

1(7710 + 1)1711 

-I- vWTTa 

m_imo(i77i + l) 5 


where we have introduced the shorthand notation m = 
(m_i,mo, TOi), and the dependence of Am on space and 
time {i,t) is implicit. We have also defined the following 
m dependent terms to simplify the writing 

Dm = -^(m-i-I-mo-I-mi)(m_i-I-mo-I-mi — 1) 

+ C /2 ^-[(mi — m_i)^ — mi — m_i]-I-mimo 

-I- m_imo^ —/r(m_i-I-mo-b mi) (13) 

S^~ = U 2 \/mi{mo + l)(mo -b 2)m_i, 

S'm''' = U 2 \/{mi + l)mo(mo - l)(m_i -b 1). 


To obtain the low energy spectrum, we follow Ref. 1^ . 
and expand the wave-function coefficients about the 
mean-field solution 


Am(*A)«e-^'>‘(AW+AW(*,t)), (14) 

where Eq is the ground state energy, Am^ is the mean field 
solution and Am^ are the fluctuations. We then expand 
the space and time dependence of Am^ into plane wave 
states to obtain 


!<(>*) = Yu Am — imomi (qt)|m_i,mo,mi) (10) 

m_i ,7710,177-1 

and the coefficients explicitly depend on the site index i 
and the time t. We now wish to variationally minimize 


AW(*,t) = (15) 

where Wk is the low energy dispersion we are trying to 
calculate. Inserting this ansatz into Eq. [1^ and keeping 
only terms linear in Uk and Uk, we obtain 
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^k'^k,m (^m -^o)'^k,m '^k,(m_i —l)(mo+2)(mi —1) *^m~^'^k,(m_i + l)(mo —2)(mi + l) (16) 

2^^k,(77i_i —l)mQmi + \/r^'l/'mo'“k,m_i(mo-l)mi + ,m_imQ(mi —1 )^ 

- t (^y/m-i + >k ,(m_i + l)momi + Vmo + ^k,m_i (mo + l)mi + VWTTiplnl* ^k,m_imo(mi + l)^ 


tlk (x/ 


TO_i + IC/^ + x/mo + IC/™ + x/jiT-i + 1C7„ / in')- 

-L m_i (m._i + l)momi '' ^ ^^t-o m_ i (mo + l)mi '' -*- m_imo(77ii + l)y 


r 


where we have introduced tpa'^ = 2z(aa), ipa'^* = 2-z(aQ) where a 7 7^ (5 £ {—1,0,1}. For example for a = 0 
[the (0) here denotes an evaluation with respect to the this would yield 
mean field ground state Am^], 

7k = 2 {cos{kx) + cos(fcy) + cos(fcz)) (17) 

and 

= J2 ^ 


(18) 


( A*i0) I 4 ( 0 ) 


= J2 1 

m' 




(19) 


{/t*( 0 ) I /l( 0 ) 


Similarly, for the vs, we obtain 


~ ^kUk rn — (-^m — -Eo)fk,m + ‘5'm ^k,(m_i-l)(mo+2)(mi-l) + ‘^'m ''-’k.(m„i + l)(mo-2)(mi + l) 


( 20 ) 


- t (x/TO-lV'™^i<,(,„_i-l)momi + x/^V^mo^m.^mo-llmi + <m_ imo(mi-1) 

- i (x/^-l + l^)^).>k,(m_ i + l)momi x/^0 “f IV’mo ^k,m_i (mo + l)mi x/^1 “1“ I'l/'mi ^k,m„imo(mi + l)) 

- t7k (x/™-l + + + VtoO + l^mo^m!_i(mo + l)mi + v'wi + llOi („i + l)) ' 


( 21 ) 


where 

= '^Vr<Tl 

m' 

X (" II' I 

^ 

m' 

V ( -I- II* ') 

^ ^k,(m'^ — l)m'^m'^ ^^(m'^-l)m'^m,gk,in'J ■ 

Because of the mean field terms, the equations of mo¬ 
tion for the Uk are coupled to those of the Wk- These 
equations can be readily solved numerically and have the 


r 


property that for every eigenvalue Wk, — Wk is also an al¬ 
lowed eigenvalue. In the following section, we discuss in 
detail the lowest collective modes Wk for repulsive and 
attractive 172 across the entire superfluid-Mott insulator 
phase diagram. 


V. EXCITATIONS 

A. Anti-ferromagnetic interactions 

We begin by discussing the case of anti-ferromagnetic 
interactions t/2 > 0. In Fig. [3] we plot the evolution of 

























FIG. 3: Evolution of the excitation spectrum from the nematic (polar) superfluid into the singlet phase. Throughout we hx 
U 2 = 0.1?7o and vary t and /r such that we access the deep superfluid (left), superfluid-Mott phase boundary (center) and the 
deep Mott phase (right). The polar superfluid has two linearly dispersing modes corresponding to density and spin modes, 
originally derived in Refs. ii As the superfluid-Mott phase boundary is approached, the sound speeds associated with the 
density and spin modes become identical. Additionally, there is a quadratically dispersing gapped, particle-hole mode and a 
non-dispersive mode. This non-dispersive mode corresponds to the breaking of a local singlet pair and has no analog in the 
spinless case. This is therefore a new feature of the spin-1 Bose gas. Deep in the Mott insulator, the singlet pair breaking mode 
has the lowest energy and the particle-hole mode is pushed to higher energies. 


the low lying excitation spectrum across the superfluid- 
singlet Mott phase diagram. Throughout, we fix U 2 = 
O.IC/q and vary t/Uo and fi/Uo with a density of n Ri 2, 
in order to access different regimes of the singlet Mott 
insulator-superfluid phase diagram. 

In the superfluid phase, we recover the three linearly 
dispersing modes; one density mode, and two degener¬ 
ate spin modes, as originally described by the authors 
of Refs. H, H for the continuum spin-1 Bose gas. The 
sound speeds associated with the density Cd and spin 
modes Cs are proportional to '/Uq and '/U 2 respectively, 
and are typically vastly different in atoms such as ^^Na 
(Ref. [i). As the Mott transition is approached however, 
density fluctuations are suppressed and consequently Cd 
approaches Cs monotonically. The sound speed associ¬ 
ated with the spin mode remains relatively unaffected, 
as it is not directly related to the compressibility. As 
discussed in Sec.IIIA, the transition from the superfluid 
to the singlet-Mott insulator is first order, and is charac¬ 
terized by an abrupt jump in the superfluid order param¬ 
eter. Thus at the Mott transition, the sound velocity also 
shows a discontinuous jump to zero (not shown). This is 
in contrast to the spinless case (or the ferromagnetic spin- 
1 gas or the nematic Mott-superfluid transition), where 
the sound speed either remains finite (at the Mott tip 
0 ) or goes to zero continuously, following the super¬ 
fluid order parameter. 

Near the Mott transition, but still in the superfluid 
phase, we additionally find two gapped modes: a quadrat¬ 
ically dispersing particle-hole mode and a non-dispersive 
mode which corresponds to the breaking of a singlet pair. 
The latter mode is absent in odd integer Mott lobes, and 
is a new feature of the spin-1 gas, with no analog in the 
spinless case. 

In Fig.m we plot the particle-hole gap Aph and the sin¬ 
glet gap As across the superfluid-Mott phase boundary. 
The particle-hole gap evolves non-monotonically as the 
superfluid-insulator transition is crossed, going to zero 
at the phase boundary, signaling a quantum phase tran¬ 
sition. For a second order quantum phase transition, the 



FIG. 4: Evolution of the singlet gap (solid) and the particle- 
hole gap (points) across the superfluid n = 2 singlet 

Mott insulator boundary for anti-ferromagnetic interactions. 
(U 2 /U 0 = 0.1). The singlet gap is completely independent of 
the hopping in the Mott insulator and and can be obtained 
by diagonalizing the on-site Hamiltonian. The particle-hole 
gap discontinuously goes to zero as the Mott-superfluid phase 
boundary is approached from the superfluid side, indicative 
of a first order transition. It approaches 0.5[/o deep in the 
Mott insulator phase. 


closing of the gap obeys the same power laws on either 
size of the transition scaling as A ^ |t — tcT where tc is 
the critical point, and u is the correlation length critical 
exponent. This is no longer true for a first order transi¬ 
tion; as shown in Fig. IH the gap closes discontinuously 
when the transition is approached from the superfluid 
side, but continuously if approached from the Mott side. 

Deep in the n = 2 Mott insulator, the gap approaches 
Aph —>■ ?7o/2 as f —0. By contrast, within the Gutzwiller 
approximation, the singlet gap is independent of the hop¬ 
ping in the Mott insulator, and can be readily estimated 
by diagonalizating the onsite part of the Hamiltonian in 
Eq. [TJ Within the Gutzwiller theory, the singlet and 
particle-hole modes do not couple to one another, and 
there is no hybridization gap as the particle-hole mode 
crosses the singlet mode. In reality, quantum fluctua¬ 
tions will couple the singlet and particle-hole modes and 
the singlet gap will depend on t. This goes beyond the 
Gutzwiller approach and has been studied by Imambekov 
et al. 0 and Snoek and Zhou , and is not reproduced 
here. Their main result is that the singlet gap indeed 














(a) 
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(b) (c) 



FIG. 5: Evolution of the excitation spectrum from the ferromagnetic superfluid into the Mott insulating phase. Throughout 
we fix U 2 = —O.IC/q and vary t and /r such that we access the deep superfluid (left), superfluid-Mott phase boundary (center) 
and the deep Mott phase (right). The ferromagnetic superfluid has one linearly dispersing mode corresponding to density 
fluctuations and a quadratically dispersing gapless and a gapped spin mode, corresponding to fluctuations of the spin in the 
direction parallel and perpendicular to the easy axis. As the superfluid-Mott phase boundary is approached, the spin gap 
associated with the spin mode goes to zero. Additionally, there is a gapped, quadratically dispersing particle-hole like mode. 
The sound speed associated with the density mode vanishes at the transition and deep in the Mott insulator, there are two non¬ 
degenerate gapped modes corresponding to particle and hole-like excitations. The gapped spin mode is effectively non-dispersive 
in the Mott insulator as its effective mass is proportional to t which is exponentially small. 



FIG. 6: Evolution of the particle-hole gap across the 
superfluid- n = 2 ferromagnetic Mott insulator boundary for 
U 2 = -O.i;7o. 

varies with hopping and jumps to zero at the singlet- 
nematic Mott transition, which is not captured within 
our theory. In general, the strong coupling theory of 
the even n Mott lobes can be described by a constrained 
quantum rotor model 

In the nematic (odd n) Mott lobes, and the ferromag¬ 
netic Mott lobe for small \U 2 \/Uq, the low energy de¬ 
scription corresponds to a ferromagnetic, spin-1 bilinear- 
biquadratic J — K spin Hamiltonian (see Eg. [22] below), 
originally derived by Imambekov et al. [lOj and Snoek 
and Zhou [1^. In the nematic Mott insulator, the low 
energy excitations are linearly dispersing nematic waves, 
whose spectrum was studied in detail by Imambekov et 
al. 0, and is not reproduced here. This theory is be¬ 
yond the naive Gutzwiller approach we develop above, as 
the Gutzwiller theory does not contain any spin fluctua¬ 
tions in the Mott lobes. Below we study this biquadratic 
spin model on the ferromagnetic side, and compute the 
spin wave spectrum and discuss its evolution as a func¬ 
tion of t. 


B. Ferromagnetic interactions 

We now turn to the case of ferromagnetic interactions 
U 2 < 0 which is shown in Fig. |SJ Deep in the superfluid 


phase, we once again recover the excitation spectrum de¬ 
rived for the spin-1 Bose gas in Refs. Si The low lying 
spectrum displays a single linearly dispersing mode cor¬ 
responding to density excitations and two quadratically 
dispersing modes corresponding to the spin excitations 
about the ferromagnetic ground state. One of the modes 
has a free-particle spectrum which corresponds to spin 
waves while the other mode is gapped, and corresponds 
to “quadrupolar” spin fluctuations 


Deep in the Mott insulator, there are two non¬ 
degenerate low lying excitations which correspond to par¬ 
ticle and hole like modes respectively. These modes are 
also present in the spinless case. At the Mott transition 
one of the gaps vanishes, signaling the phase transition 
into a superfluid and the other gap remains finite across 
the phase boundary as shown in Fig. [51 Near the phase 
boundary but on the superfluid side, the effective mass 
of the spin mode decreases, as in a lattice, the effective 
mass is proportional to the hopping, which scales expo¬ 
nentially with the lattice depth. Furthermore, the spin 
gap associated with the quadrupolar spin mode vanishes 
at the transition. At the ferromagnetic transition the su¬ 
perfluid density vanishes, leading to a vanishing phonon 
velocity. 


Note that unlike the spinless gas however, the ferro¬ 
magnetic Mott insulator has long range spin order, and 
thus according to Goldstone’s theorem, has a gapless 
mode corresponding to spin excitations. However this 
mode is not captured within the naive Gutzwiller ap¬ 
proach as spin fluctuations in this theory are tied to the 
condensate order parameter, and therefore vanish at the 
Mott transition. To capture spin fluctuations in the insu¬ 
lating phase therefore, we augment the Gutzwiller theory 
with a Schwinger boson mean-field theory for the spin, 
which is described next. 
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1. Schwinger boson mean-field theory 

In the following subsection we determine the low lying 
excitations in the ferromagnetic Mott phase using the 
Schwinger boson mean field theory (SBMFT). We 
will focus on the n = 1 Mott lobe for simplicity, but 
our results are also relevant to any value of n with a 
change of the spin model parameters. Working in the 
limit Uq ^ 1 1 / 2 1 t, we can apply standard perturbation 
theory in t/{Uo + 51 / 2 ) (where g is a small integer) which 
yields the spin-1 biquadratic model [T^. 

Hjk = Y. ) . ( 22 ) 

The parameters J > 0 and Kj^ > 0 are related to 
t,C/o,t /2 via Eq. 22 of Ref. [13. Focusing on the ferro¬ 
magnetic ground state we can assume that the K term is 
not sufficient to destroy the long range ferromagnetic or¬ 
der. Thus, we can treat the biquadratic term at the mean 
field level, which yields (S^ • Sj)^ —>• —d’fj + 2<I)ijSi • S^-, 
where = (Si-Sj), and therefore this decoupling essen¬ 
tially just renormalized the nearest neighbor ferromag¬ 
netic interaction to J -I- K^ij. We generalize the spin 
symmetry from SU(2) to SU(IV) and introduce Schwinger 
bosons through = b\^Smnbim where Smn are the gen¬ 
erators of SU(7V) and m,n = 1,...,7V. The b opera¬ 
tors must satisfy the constraint = NS, and 

the SBMFT is exact in the limit TV —>■ 00 . We stress 
that the bosonic spinon operators bin are not the bosonic 
operators in the Bose-Hubbard model in Eq. |T] Follow¬ 
ing standard SBMFT techniques we decouple the spin- 
spin interaction in the ferromagnetic channel through 
S, . S, =: : /TV - 5^, where T,, = Enblb.n, 

and : • • • : denotes normal ordering. Solving the SBMFT 
equations with TV = 2 and S' = 1 at zero temperature 
yields a ferromagnetic ground state [5l| with a bosonic 
excitation spectrum given by 

Wk = z{J + K) - {J + K)ji, (23) 

(J + iV)|k|2. (24) 


We recover the expected quadratically dispersing fer¬ 
romagnetic spin waves with an effective mass m* = 
1/(2[J -I- K]) ^ (in units of Ti = 1). There¬ 

fore, we conclude that the correct excitation spectrum 
in the ferromagnetic Mott lobes have gapless Goldstone 
modes which disperse quadratically and cannot be cap¬ 
tured within the Gutzwiller approach. We emphasize 
that Eq. [22] is only valid deep in the Mott insulator 
and not near the transition where the truncation of basis 
states needed to arrive at this equation is no longer valid 
due to the vanishing of the particle-hole gap. 


VI. EXPERIMENTAL IMPLICATIONS 

The continuum physics of the spin-1 Bose gas has been 
well studied experimentally, and the phase diagram is 
well established [3,1 Is^ - issjl . By contrast, the phase di¬ 
agram of the lattice spin-1 gas has received relatively 
little attention from the experimental community, de¬ 
spite the plethora of interesting phases and phase transi¬ 
tions present in this model. As we have shown here, the 
strongly correlated spin-1 superfluid and Mott regimes 
have many distinct features that are absent in the well 
studied spinless Bose Hubbard model Indeed it will 
be extremely interesting to study the evolution of the 
nematic and singlet order parameters in a strongly cor¬ 
related spin-1 gas with anti-ferromagnetic interactions, 
as in ^^Na. Importantly, the evolution of these order 
parameters reveals a first order quantum phase transi¬ 
tion near the superfluid-singlet Mott insulator, which has 
been confirmed numerically in ID. It will be very exciting 
if the predicted first order nature of the quantum phase 
transition can be explored experimentally. 

The excitations in the spin-1 gas are also strikingly 
different from their spin-0 counterparts. For example, 
unlike the spinless Mott insulator, which is truly feature¬ 
less, the ferromagnetic Mott insulator has quadratically 
dispersing spin waves, corresponding to long range spin 
order. The excitations in the weakly interacting super¬ 
fluid limit of the spin-1 ferromagnetic gas were recently 
explored by Marti et al. , where a spin wave was ex¬ 
ternally imprinted on to the cloud and its coherent evolu¬ 
tion was subsequently imaged. This method can also be 
applied in the Mott insulating regime to explore the spin 
wave spectrum in the ferromagnetic Mott lobe. Particle- 
hole like excitations couple to density fluctuations which 
are readily generated using modulation spectroscopy 
or Bragg spectroscopy |58j |. For anti-ferromagnetic in¬ 
teractions, the low energy excitations are nematic waves, 
which are linearly dispersing in the superfluid and ne¬ 
matic Mott insulator. As is known from the theory of 
liquid crystals, the nematic tensor couples to the polar¬ 
ization of the incoming light beams, leading to optical 
birefringence, which can be used to probe nematic order 
and nematic waves (bfij . 

VII. CONCLUSIONS AND OUTLOOK 

To conclude, in this paper we have presented a compre¬ 
hensive mean-field description of the static and dynamic 
properties of the superfluid-Mott insulator transition in a 
spin-1 Bose gas. A key distinction in our work from pre¬ 
vious works on this subject is our focus on the evolution 
of the important order parameters for the ferromagnetic 
and anti-ferromagnetic interactions, namely the spin in 
the ferromagnetic case, and the singlet and nematic order 
parameter in the anti-ferromagnetic case. 

Additionally, we have described the low lying excita¬ 
tion spectrum of the spin-1 gas across the entire phase di- 
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agram for positive and negative C/ 2 - For C /2 > 0, we have 
discussed the evolution of the singlet and the particle- 
hole gap, which can be probed using modulation spec¬ 
troscopy [13 • The singlet gap is a new feature of the 
spin-1 gas, and has no analog in the spinless case. We 
have discussed the limitations of the Gutzwiller approach 
in that it neglects the quantum fluctuations which cou¬ 
ple the singlet and particle-hole modes. This would make 
the singlet gap vary as a function of tunneling, eventually 
going to zero at the singlet-superfluid transition. 

On the ferromagnetic side, we have discussed the evo¬ 
lution of the quadrupolar spin gap and the particle-hole 
like excitations across the phase diagram. Furthermore, 
we have pointed out a shortcoming of the Gutzwiller ap¬ 
proach in describing spin fluctuations in the Mott in¬ 
sulator. Unlike the spinless Mott insulator, the ferro¬ 
magnetic Mott insulator is not featureless but rather is 
characterized by long range spin order. However, within 
this theory, spin fluctuations are tied to the condensate 
order parameter. Therefore, this theory accurately cap¬ 
tures the spin modes in the condensate and reproduces 
the Bogoliubov spectrum at weak coupling. However, in 
the Mott insulator, where the condensate order param¬ 
eter vanishes, spin fluctuations are frozen out, and as a 
result, the ferromagnetic Mott insulator does not have 
any spin fluctuations, in violation of Goldstone’s theo¬ 
rem. To overcome this limitation, we have presented a 
Schwinger boson mean-field theory, which retains spin 
fluctuations in the Mott insulator and yields an addi¬ 
tional free-particle like mode with an effective mass m* 
which varies like Uq within the Mott insulating phase. 
This mode is a new feature of the spin-1 gas, and can be 
probed using magnon interferometry [s^. 

Theoretically, the Gutzwiller approach developed here 


serves as a natural starting point for exploring more com¬ 
plicated Hamiltonians where the single particle physics 
involves a coupling between spin and kinetic degrees of 
freedom, such as the spin-orbit coupled Bose Hubbard 
model. The interplay between large spin and spin-orbit 
coupling can lead to simultaneous nematic, ferromagnetic 
orders with broken translational symmetry (s^ or exotic 
spin models with novel ground states even at the classi¬ 
cal level [Hiiii. In the presence of single particle degen¬ 
eracies such as those introduced by spin-orbit coupling, 
or artificial gauge fields, the absence of a bosonic Pauli 
principle severely limits exact numerical approaches, and 
only small system sizes can be accurately simulated nu¬ 
merically. Extending the Gutzwiller method to study 
the mean-field physics of these large spin, spin-orbit cou¬ 
pled models is therefore imperative [621 l63l| , and serves 
as a useful starting point for exploring the role of quan¬ 
tum fluctuations, the breakdown of mean-field theory and 
other strongly correlated effects, such as the fermioniza- 
tion of bosons in flat bands . Importantly, this mean- 
field theory can be systematically generalized to incor¬ 
porate fluctuation effects by solving the system exactly 
for small clusters, coupled by mean fields or by supple¬ 
menting the Gutzwiller method with Schwinger bosons as 
done here to correctly capture low energy spin physics. 
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